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Abstract 

To obtain high-quality results from lattice QCD, it is important to use operators that produce 
good signals for the quantities of interest. Free-form smearing is a powerful tool that helps to 
accomplish that goal. The present work introduces a new implementation of free-form smearing 
that maintains its usefulness and reduces its computational time dramatically. Applications to 
the mass spectrum of B, Bg, Be and bottomonium mesons show the effectiveness of the method. 
Results are compared with other lattice QCD studies and with experimental data where available. 
The present work includes the first lattice QCD exploration for some of these mesons. 
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I. INTRODUCTION 


In the usual Lagrangian formulation of lattice quantum chromodynamics (QCD) physical 
quantities are extracted from Euclidean-time correlation functions of operators which serve 
as interpolating helds for the states of interest. In general, many states contribute to a given 
correlation function. Contemporary applications of lattice QCD carried out to study radial 
and orbital excitations of hadrons can beneht from methods that can enhance signals in a 
correlation function from particular states selectively. This is what free-form smearing aims 
to do. 

A common approach to operator building is to begin with a gauge-invariant operator 
dehned at a single lattice site (or a small number of sites) that has the desired quantum 
numbers, and then to smear that operator spatially through an iterative Gaussian function 


0 


ll-l3|. At the time this method was developed the goal was to suppress excited-state contri¬ 
butions to the correlation functions. As pointed out in Ref. jd], this approach can require 
many iterations, especially on fine lattices, and the resulting operator is typically not Gaus¬ 
sian shaped. This led the authors of Ref. ^ to propose a method called free-form smearing, 
which can produce a gauge-invariant operator of any desired shape. The new method still 
uses an iterative Gaussian approach but combines it with a reweighting formula involving 
a function that the user can choose arbitrarily. This allows one to enhance or reduce the 
contribution of particular states to the correlation function. Analyzing simultaneously cor¬ 
relators constructed from operators with different free-form smearing profiles allows for a 
more reliable separation of contributions of different states. At present, free-form smearing 
is practical at the source and not at the sink, but it still provides a clear improvement over 
traditional Gaussian smearing, as will be demonstrated in the present work. 

Our implementation of free-form smearing differs from that of Ref. jd] by avoiding the 
iterative Gaussian step entirely. On an x T lattice, this provides a speedup of order L. 
This version of free-form smearing is introduced in Sec. UTl 

For a demonstration of the power of free-form smearing for physical quantities, the re¬ 
mainder of this paper presents a study of meson masses containing at least one bottom 
quark. Section HTTI explains how free-form smeared operators are built in lattice nonrelativis- 
tic QGD (NRQGD), using hydrogen-like wave-function shapes to target radial excitations 
and orbital excitations in addition to the ground state. Section [IV] provides details about 
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the ensemble of configurations, about the NRQCD action for bottom quarks, and about the 
relativistic actions used for up/down, strange and charm quarks. 

The masses and mass splittings obtained for the bottomonium spectrum are presented 
in Sec. El including some states that have not yet been observed experimentally nor in 
previous lattice studies. The B, Bg and B^. mesons are discussed in Sec. IVTl where the 
Be spectrum in particular provides predictions for upcoming experiments. Previous lattice 
studies have considered some systematic effects that cannot be addressed in the present work 
and therefore provide a valuable context, so we make direct comparisons where possible. In 

rum can be found in Refs, and 


particular, lattice studies of the bottomonium spec 
lattice studies of bottom mesons appear in Refs. [l4 




II. FREE-FORM SMEARING 


Quark smearing is a technique used frequently in lattice QCD spectroscopy to improve 
the signal obtained from correlation functions. The general idea is to smear the quark helds 
until they form an operator that resembles the wave function of the desired bound state. 

nn 

One popular method [1H3[ is to smear the quark field into a Gaussian shape by iteratively 
applying the discrete gauge-covariant Laplacian operator A as given by 


ij{x) = (l + , (1) 

where n is the number of iterations and a is the smearing parameter. Note that the field 
is only smeared in the spatial directions and not the temporal direction. This approach 
is straightforward to implement, improves the ground-state signal, and suppresses excited 
states. 

Another smearing method that is commonly used in lattice NRQCD is to £x the gauge 
links to Coulomb gauge and explicitly give the quark held an arbitrary shape using the 
formula 


y 


( 2 ) 


One can tune the function /(x 


y) to obtain optimal signals not only for the ground state. 
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but for excited states as well. By choosing a shape that closely matches the wave function 
of an excited state, the ground state can be suppressed and a cleaner excited-state signal 
obtained. The convolution can be efficiently implemented by a fast Fourier transform, and 
gauge hxing is required because Eq. ([2]) is not gauge invariant. 

Free-form smearing jd] combines the advantages of the previous methods. It allows the 
quark held to be smeared to an arbitrary shape while retaining gauge invariance, without 
the need for gauge hxing. Free-form smearing was initially applied to relativistic quarks, 
but we apply it for the hrst time to nonrelativistic heavy quarks. In our initial study of 
free-form smearing applied to bottomonium, a reduction in statistical errors relative to the 
gauge-hxed method was observed |^. 

In the original implementation of free-form smearing jd], the quark held ^ single 

lattice site x was Gaussian smeared as in Eq. ([T]) so that gauge links reach from site x to 
all other spatial sites on a given time slice. The free-form smeared held is given by the 
reweighting formula 


= ^f{x-y) 


'4’x{y) 


^x{y) 


( 3 ) 


where f{x — y) is an arbitrary function, i^xiy) is the component of a Gaussian smeared held 


which gauge transforms at x with the quark held ^|J located at ?/, and 


'^x{y) 


is the 


ensemble average of the norm of the Gaussian smeared held. The norm is dehned by 


iJxiy) 



{^l{y)^x{y)^ 


1 


( 4 ) 


where the trace is over spin and color indices. Note that the Gaussian smeared held 'ijj{x) 
in Eq. ([T]) is obtained from iJx{y) by summing over y. Any number of free-form smeared 
helds 'ip{x) can be generated by reusing the same held iJx{y) and choosing a diherent shape 
f{x-y). 

Free-form smearing is fairly insensitive to the Gaussian smearing parameters a and n. 
The parameter n should be chosen large enough so that gauge links reach every spatial lattice 
site, and a must be chosen so that f < |, otherwise Gaussian smearing breaks down. For 
a lattice with L sites in each of the spatial directions, n = Gaussian smearing steps are 
required to reach all sites. This large number of Gaussian smearing steps is a computational 
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FIG. 1: Illustration of the minimal-path method of free-form smearing in two dimensions. The 
first, second and third iterations are shown from left to right, respectively. The source position is 
in the center of each sketch and links are multiplied outwards iteratively, forming a shell. Different 
link paths that come to the same site are summed; the number of minimal paths for each site is 
indicated. 


burden for the original implementation of free-form smearing. 

A much more computationally efficient alternative to the Gaussian version of free-form 
smearing is to replace the initial step ([T]) by minimal paths, that is, the shortest paths 
along links from one lattice site to all other spatial sites. Recall that Gaussian smearing 
generates many link paths from one site to all other spatial sites, but that the links are 
also multiplied by a factor a/n which must be less than 1/6. Therefore, since the Gaussian 
method multiplies the links by a small factor, the shortest paths tend to dominate. 

The new free-form link algorithm starts at a site x and multiplies links outward to all 
nearest neighbors, and then to all next-nearest outward neighbors, as illustrated in Fig. [T] 
Paths along different links that lead to the same site are added together, thus resulting in 
a sum of minimal paths. This continues until all spatial sites have been reached and every 
link has been multiplied exactly once. The result is all minimal paths from the original site 
X to all other spatial sites ?/, 


Uv) = Z U{x ^ y)'il:{y) . (5) 

shortest paths 

Free-form smearing can then be implemented as in Eq. (I3]). Whereas the Gaussian version of 
free-form smearing requires 0{L‘^) link multiplications, the minimal-path method uses 3L^ 
link multiplications, which is exactly the number of spatial links on a given time slice. 

Link smearing can also be used prior to the Gaussian smearing or minimal-path step. 
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FIG. 2: Comparison of effective masses for the Gaussian and minimal-path versions of free-form 
smearing, using the same function f(x — y), applied to a bottomonium pseudoscalar correlation 
function. The bottom panel shows a close-up view of the ground-state plateau. The minimal-path 
method is much more computationally efficient. 


In the present work, stout-link smearing is used with parameters p = 0.15 and Up = 10 as 


dehned in Ref. 


22| . Therefore, our minimal-path method actually uses the shortest path of 


stout links, which effectively contains longer link paths as well. Link smearing is useful in 
reducing some excited-state contamination, although the most signihcant improvements for 
the present work are made by careful selection of the smearing shape f{x — y). 

Figure |2] shows a direct comparison of the Gaussian and minimal-path free-form smearing 
methods for the effective mass of a bottomonium S-wave correlation function. The same 
wave-function shape f{x — y) and stout-link smearings are used in both cases, and the 
Gaussian smearing parameters are n = 64 and « = 0.15 from Eq. O- There is virtually no 
difference between the results of the two methods shown in Fig. [2] and the same is true for 
other operators as well; what slight differences can be found (e.g. in a P-wave correlation 
function) are only an issue of hne-tuning in f{x — y). 

For the present study, smearing the held over the whole spatial volume is appropriate 
since the volume is not very much larger than the hadron size. This also conforms to the 
way Gaussian smearing, wave-function smearing and the original free-form smearing have 
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been done in the past [l-^. 


III. LATTICE OPERATORS 


For this study, smearing is applied to the bottom quark, leaving the antiquark unsmeared. 
This is expedient because the bottom-quark propagator is obtained from lattice NRQCD and 
is less computationally expensive than the relativistic propagators for the up/down, strange 
and charm quarks. For bottomonium, the resulting wave function places an anti-bottom 
quark at the origin surrounded by a bottom quark. For bottom mesons, an anti-up/down, 
anti-strange or anti-charm quark is surrounded by a bottom quark. A physically intuitive 
picture would instead place the center of mass at the origin but, since we take a zero- 
momentum projection of the meson, the location of the center of mass is irrelevant; only the 
distance between the quark and antiquark matters. 

Hydrogen-like (i.e. Coulomb potential) wave-function sha^s have been used successfully 
within the gauge-fixed smearing method in lattice NRQCD p and we use them here with 
free-form smearing as well. Different shapes are used depending on the intended orbital 
angular momentum, and nodes are included to optimize the operator for radial excitations. 
Here is a list of the basic smearing shapes used in this paper: 


e “0 


S-wave: f{x — |/) = < e “o (r — 6) 

__ r 

e “0 (r — c)(r — b) 


P-wave: fi{x - y) = 
D-wave: fij{x-y) = 


“0 Xi 


e “0 Xi (r — h) 

e ^ {xiXj — ^dij{xl x\ + 

e“'^ {xiXj - yij{xl + xl + xl)){r -b) 


F-wave: fijk{x - y) = XiXjXk e “o 

_ ^ 

G-wave: fijki{x - y) = XiXjXkXi e “o , 


( 6 ) 

(7) 

( 8 ) 

(9) 

( 10 ) 


where Xj is defined to be periodic, but r is simply defined as the shortest distance between 
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sites X and ?/ in a periodic box: 



2 

min ’ 


( 12 ) 


( 11 ) 


The radius and nodal parameters (oq, c) are tuned to optimize the signal of the ground- 
state or radial excitations. Free-form smeared operators are built according to 


x{x)ii(x) = x(x) fi/(i - y )— 



(13) 


where ijj is the quark, x is the antiquark and Qf{x — y) is the free-form operator. Table [I] lists 


all operators used in this paper by their lattice irreducible representations and continuum 
quantum numbers. It provides a thorough coverage of quantum numbers for S, P and D 
waves plus hve operators that offer the simplest exploration of some of the F and G waves. 


Ill |. additional F and G operators would duplicate some of the 


As explained in Ref. 


quantum numbers that already exist in this table, so they are omitted here. Sink operators 
are unsmeared and are given by the covariant derivative operators listed in Tables II and 
III of Ref. [n|. 

The optimized free-form smearing parameters corresponding to Eqs. used in this 

study are given in Table UTl The parameters (oq, b, c) are tuned by hand to optimize the signal 
in a correlation function for the ground state, the hrst excited state and, for bottomonium S 
waves, the second excited state. These optimizations improve the accuracy and precision of 
the spectrum. The effective mass plots shown in Figs. [3] and |4] demonstrate the effectiveness 
of free-form smearing to obtain clean signals for ground states and excited states. Figure [3] 
shows that the plateau begins at small Euclidean times for a bottomonium D-wave ground 
state and for the hrst radial excitation. In particular, note that the effective mass for the 
hrst radial excitation shows no contamination from the ground state. Similarly, Fig. |4]shows 
small-time plateaus for meson S-wave, P-wave and D-wave ground states. 

Free-form smearing shrinks the contamination from unwanted states, thereby emphasizing 
the one state of interest, but it does not eliminate the contamination entirely. Diherent 
choices for the smearing parameters can be used to emphasize each state separately, and 











TABLE I: The factor that defines the quantum numbers of the operator in Eq. (|13l) and incorporates 
the free-form smearing functions from Eqs. dSD-dini). Column 2 shows only the J value that is 
expected to dominate in each case; for other possible J values see Table I of Ref. [ll|. 
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{(/ll — / 22 )/V^, (/ll + /22 — 2 / 33 )/\/ 6 } 

{/23, /3I) /12} 

+ /l20'2 + flsO'S, f2lO'2 + /220'2 + /230'3, /si^"! + /320‘2 + fssO's} 
{(/230-1 - /i 30-2)/\/2, (/23Cri + /3lO-2 “ 2/i2a3)/\/6} 

{(/22 — + /l3Cr3 — fl20'2, {fss “ /ll)cr2 + f2lO'l “ /23Cr3, 

(/ll — /22 )o‘ 3 + /320‘2 — /siCTl} 

13^1 + f31<^2 

2 /i 20'2 ~ 2/13(73, 3f220'2 — 2/23*73 ~ 2/21(71, 


22)<73 + 2 / 31(71 — 2 / 32 ( 72 } 

/l33) /233 “ /2II, /3II “ /322} 

/33i)<7i + (/332 “ /ll2)<72 + (/ll3 ~ f223)<^3 
' /3332, / 333 I — / 1 II 3 , / 1 II 2 — / 2221 } 

/3332)<7i + (/ 333 I — /i 113)<72 + (/lll2 “ /222l)<73 


then a simultaneous fit to several correlation functions will be stable and reliable. For the 
present work, nonoptimized correlation functions are also included in the multistate fits 
and they further stabilize the results. For bottomonium, the nonoptimized correlators are 
calculated using unsmeared operators as given in Tables II and III of ll|. For B, Bg and 
Be ground-state profiles, small radial parameters oq = 0.5 or oq = 1.0 are used as the 
nonoptimized choices. 

Chronologically, the calculation for Be was started before we had developed the minimal- 
path method, so it uses the Gaussian version of free-form smearing with parameters n = 64 
and a = 0.15 from Eq. ([T]). The minimal-path method of free-form smearing is used for all 
B, Bg and bottomonium calculations. Throughout this work, stout links with parameters 
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TABLE II: Free-form smearing parameters for B, Bg, Bf. and bottomonium. Parameters were 
optimized by trial and error. In addition, nonoptimal parameters (ao = 0.5 and 1.0) are used for 
B, Bg and B^ as well as unsmeared operators for bottomonium, as discussed in the text. 


bottomonium 


Be 



Bg 



B 


L 

ao 

b c 

L 

ao 

b 

L 

ao 

b 

L 

ao 

b 

S 

1.4 


S 

0.5 


S 

0.5 


S 

0.5 



1.6 



2.0 



4.5 



5.0 



2.6 

2.6 


2.2 



5.8 

5.8 


5.5 



2.8 

2.8 


4.0 

4.0 


6.0 

6.0 


6.2 

6.2 


3.0 

2.13 6.0 

P 

1.0 


P 

1.0 


P 

0.5 


P 

2.0 



3.0 



6.0 



7.0 



3.0 

4.6 


4.0 

6.0 


7.0 





D 

2.5 


D 

1.0 









3.5 

6.5 


5.0 









F 3.0 
4.0 
G 4.0 
5.0 


p = 0.15 and Up = 10 were employed. 


IV. SIMULATION DETAILS 


An Nf = 2 + 1 gange field ensemble from the PACS-CS Collaboration 
this study. The action contains an Iwasaki gauge term and clover-improved 
The lattice dimensions are 32^ x 64 and the lattice spacing 


23| was used for 


Wilson fermions. 


opAcs-cs = 0.0907(13) fm 


(14) 


was determined by the PACS-CS Collaboration 23(] using the experimental vr, K and U 
masses. The pion mass is near physical at = 156(7) MeV. At (mA')orig = 554(8) MeV, 
the kaon mass is larger than its physical value so, to account for this, a partially quenched 
strange quark is used in this study as was done in Ref. 2^ where the valence strange quark 
was retuned to match the physical 0 mass, giving a kaon mass of rriK = 504(7) MeV. 

The charm quark parameters are taken from Ref. 2^ where they were tuned using the 
Fermilab interpretation. The parameters for the relativistic quarks are given in Table IIIII 
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FIG. 3: Effective mass plots for the free-form smeared ^D 2 {E representation) bottomonium 
correlation functions tuned to optimize the ground state and the first excited state. Energy val¬ 
ues extracted from fits to the correlation functions are shown as solid horizontal lines. Dashed 
horizontal lines are the statistical uncertainties. 


TABLE III: The hopping parameters and clover coefficients for relativistic sea and valence quarks. 


Quark flavor 

K 

Csw 

sea quarks 



up/down 

0.13781 

1.715 

strange 

0.13640 

1.715 

valence quarks 



up/down 

0.13781 

1.715 

strange 

0.13666 

1.715 

charm 

0.12686 

1.64978 


The propagators for the relativistic valence quarks {u/d, s, c) are calculated using the sap_gcr 


solver from the DD-HMC code made available by Martin Liischer 


3-27|. 


The bottom quark is implemented using lattice NRQCD, including terms up to 0(v‘^ 


in the bottom-quark velocity, which corresponds to q = 1 
i > 7 according to the notation of the Appendix in Ref. 


or 1 < z < 6 and q = 0 for 
281. Following Ref. [ll|, the 


bare mass of the bottom quark [Mf, = 1.95) is taken from fitting the kinetic mass of the 
T to its experimental value, the tadpole factor is set to the average link in Landau gauge 
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FIG. 4: Effective mass plots for the free-form smeared 0“, 1'*' and 3“ {T 2 representation) Be 
meson correlation functions tuned to optimize the S-wave, P-wave and D-wave ground states. The 
P-wave operator shown corresponds to I"*'’*' when applied to quarkonium; see Sec. IVIl for discussion 
on this operator when applied to heavy-light mesons. Energy values extracted from fits to the 
correlation functions are shown as solid horizontal lines. Dashed horizontal lines are the statistical 
uncertainties. 


{ul = 0.8463), and the stability parameter is chosen to be n = 4. 

The calculation of correlation functions for bottom mesons with a nonrelativistic bottom 
quark requires some additional care with regard to the contraction of spin indices. Our 
choice for Dirac 7 matrices is 


7fc = 




74 




(15) 


which allows the Dirac indices of the relativistic propagator and the Pauli indices of the 
nonrelativistic propagator to be contracted in a simple manner. For forward time propaga¬ 
tion the top two Dirac indices are used, while for backward time propagation the bottom 
two Dirac indices are used. 

Simulating multiple sources on a given configuration can reduce statistical errors. A 
random U{1) wall source imitates multiple sources at different spatial sites on a given time 
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slice without the need to calculate many propagators. While it is easily implemented using 
the conventional smearing methods given by Eqs. ([T]) and ([2]), a free-form smeared random 
wall source is computationally expensive because every site x must be smeared independently 
and summed according to the formula 

_ N 

= (16) 
i 

where 0 < 9w{xi) < 27r is a random complex phase at each lattice site. An efficient approach 
is to calculate a partial random wall source, where only a relatively small number of spatial 
sites X is used. For this study the partial wall source has = 4^ evenly spaced free-form 
smeared sites. With the original implementation of free-form smearing, the cost of free-form 
smearing this partial wall source was comparable to the cost of calculating one NRQCD 
propagator over the full time extent of the lattice. The new implementation of free-form 
smearing reduced the cost by a factor near 100. Including more sites in the partial wall source 
did not improve the statistical errors for the case of bottomonium . Statistical errors are 
further reduced by using partial wall sources for multiple time slices. For bottomonium a 
partial wall source is built for each of the 64 time slices, and for the B, Bg and Be mesons 
a partial wall source is built for 32 time slices. To further increase the statistics, NRQCD 
propagators are calculated forward and backward in time. The correlation functions are 
binned over all time sources and over forward and backward propagation. 

The simulation energies are extracted by a simultaneous multiexponential ht to multiple 
free-form smeared correlation functions of the form 

= , (17) 

n 

where En is the fit parameter for the energy of the nth state. Each channel from Table 
[T] is fit separately. The fits are done to correlation functions with free-form smeared source 
operators as described in Eqs. l)6 |) - ffTUD and Tables [H and HIl For bottomonium, unsmeared 
source operators are used as well. All hts exclude the source time step at t = 0, and are 
typically truncated at t = 15 or t = 23 because the signal is lost in noise for larger Euclidean 
times. 
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V. BOTTOMONIUM SPECTRUM 


For bottomonium, the S-wave correlator is fitted to six exponentials, where the lightest 
three energies are identified as the ground state, first radial excitation and second radial 
excitation. The other exponentials are collectively able to account for higher excitations but 
we do not interpret them individually as the physical states. 

For the P waves and D waves, separate fits are performed to extract the ground states 
and hrst excited states. The fits for the P-wave and D-wave excited states are done using 
three correlation functions: an unsmeared operator, a free-form smeared ground state, and 
a free-form smeared excited state. The ground-state P waves and D waves are ht to the first 
two correlation functions with the excited-state optimized correlator excluded. Removing 
this correlation function from the £t significantly reduced the statistical uncertainty of the 
ground states. These more precise ground-state ht values are statistically compatible with 
the noisier ground-state ht values that are obtained by including the excited-state correlation 
function. The excited-state hts use hve exponentials (six for ^Di and ^D^Ti) and the ground- 
state hts use four (hve for ^Di and The D-wave operators ^Di and have 

a noticeable mixing with the ^Si ground state and an additional exponential is included 
in the ht for these two channels. The smallest energy is consistent with the ^Si ground 
state and all higher energy levels are assumed to be D waves. A conhdent interpretation 
would require a cross-correlation matrix between the ^Si, ^Di and ^D^Ti operators, but this 
is beyond the scope of the present study. We note that the results for ^Si also assume no 
D-wave contamination. 

The F-wave bottomonium state is ht with four exponentials to three correlation 
functions: unsmeared, free-form smeared with oq = 3.0, and free-form smeared with oq = 
4.0. To obtain a reliable result, is ht with four exponentials to the unsmeared and 

Oq = 3.0 free-form smeared correlation functions, and ^F 3 T 2 is ht with three exponentials 
to the Oo = 3.0 and oq = 4.0 free-form smeared correlation functions. The bottomonium 
G-wave ground states are obtained from four exponential hts to three correlations functions: 
unsmeared, free-form smeared with oq = 4.0 and free-form smeared with oq = 5.0. Even 
for these higher orbital angular momenta, the hydrogen-like smearing prohles from Eqs. ([9]) 
and flTOj) produced signihcantly better ground-state signals than Gaussian smearing. 

Since the bottom-quark mass is hxed by tuning the kinetic mass of the T to its physical 
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TABLE IV: Bottomonium mass splittings for IP, 25, 1^P2 and 2P with respect to 15, where a 
bar represents the spin average. Results using two different lattice spacing definitions are shown: 
Eq. (|14l) and Eq. (|19l) . The value for ^D 2 is the dimensional average of the E and T 2 lattice 
representations. 




m — niYg [Me V / 



lattice using opAcs-cs 

lattice using ajp jg 

experiment [3h| 

Tp 

437(6) 

455.0(9) 

455.0(9) 


547(10) 

569(6) 

572.5(1.3) 

I^D2 

688 (10) 

715(3) 

719.0(1.7) 

2P 

743(15) 

773(11) 

815.4(9) 


value 



the absolute masses of the bottomonium spectrum are calculated from 


he 


mc^ = 


(18) 


Using the lattice spacing scale given in Ref. 


23|. the masses of the 25, IP, 2P and ID bot- 


tomonia are all systematically smaller than the experimental values by a significant amount, 
as shown in Table IIVI However, ratios of differences having the 15 mass subtracted from 
the 25, IP and ID masses agree with experiment, which suggests using the bottomonium 
spectrum to set the scale. A new lattice spacing is defined using the spin-averaged IP — 15 
mass splitting 


a _- he ^is)sirr 

“lP-15 — no n ox 


= 0.0872(3)fm , 


(19) 


exp 


which is 4.0% smaller than the PACS-CS value from Eq. flT4|) . Reference 2^ noted that 
other methods find a lattice spacing that is up to 4.4% larger than the PACS-CS value for 
this ensemble. 


The lattice spacing reported in Ref. 23| was obtained from the light quark hadron spec¬ 
trum, which is less relevant for the case of bottomonium. Therefore, we use Eq. (IT^ to set 
the scale for the bottomonium spectrum. Even after using the scale from Eq. ([19]), the 2P 
masses remain systematically smaller than experiment. This cannot be due to contamination 
from higher excited states because that would cause the 2P masses to be larger, not smaller. 
Applying radiative corrections or nonperturbative tunings to the NRQCD coefficients q and 
including higher-order terms in the bottomonium velocity are possible ways to remove this 
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FIG. 5: Mass spectrum of bottomonium. Red bands are experimental values. Black points with 
errors bars are lattice data with statistical errors only. Grey bands are the statistical error and 
lattice spacing uncertainty, added in quadrature. Numerical values are given in Table IVTIl 
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discrepancy. In principle, the bottom-quark mass should also be retuned using this new 
lattice spacing. At present, we simply note that setting the scale with physics relevant to 
the bottom quark increases the accuracy of the bottomonium spectrum. 

The entire bottomonium spectrum below the BB threshold [with the exception of the 3P 
states, where the experimental value of Xfei(3-P) is just below the BB threshold 29|] is shown 
in Fig. |5l as extracted using chi-squared hts of free-form smeared correlation functions. For 
spin-2 and spin-3 states, where results for more than one lattice irreducible representation 
A were calculated, the dimensional average of the simulation energies. 


^^ Ea dim(A)PA 

EAdim(A) 

is our reported value for the mass. The grey bands show a combination of the statistical 
errors and a 4.0% systematic uncertainty in the lattice spacing, which come from discrep¬ 
ancies in the determination of the lattice spacing. Our work provides the hrst lattice result 
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for the bottomonium D-wave radial excitations in all channels. 

Precise spin splittings were obtained for IS", 2S, IP, 2P and ID bottomonia, shown in 
Table El The spin splittings in Table E agree well with experiment, except for the Xboi^P) 
which is larger. This discrepancy with the BABAR results 3l| is more evident in the 
spin-dependent splittings 


34.5(9) MeV 

this work 

( 21 ) 

46.0(1.9) MeV 

BABAR 

157(4) MeV 

this work 

( 22 ) 

160.0(2.2) MeV 

BABAR 


~2xf,o(l-P) + 3x6i(1-P) — Xb2(l-P) 
-2xbo{^P) - Sxbi{lP) + 5x62(1T’) 


where, at tree level, the former is proportional to the NRQCD parameter c| while the latter 
is proportional to C 3 . For Xb 2 {^P) in Eqs- (El]) and (|2^ the dimensional average of the E 
and T 2 lattice irreducible representations is used. Since our simulations only use tree-level 
coefficients where C 3 = C 4 = 1, the P-wave spin structure could be improved by tuning C 4 
until Eq. fET]) agrees with experiment. Reference {q] noted that increasing C 4 > 1 had the 


effect of decreasing the Xboi^P) mass relative to the spin-averaged mass while doing 
little else. 

The D-wave splittings have not yet been experimentally observed, but our results in Ta¬ 
ble Elare smaller than the predictions given in Ref. 10 1, also shown in Table E Reference Q 
proposed a method to reduce systematic effects from C 3 and C 4 dependence. When applied 
to our data, this method produces results that are consistent within statistical uncertainties 
with the simple difference shown in Table E 


VI. B, Bs AND Pc SPECTRUM 

In contrast to the case of bottomonium, charge conjugation is not a helpful quantum 
number for bottom mesons. Whereas the I"*"” and operators in Table [I] couple to 


separate quarkonium s 


light meson states 32 


ates, those same two operators each couple to a mixture of heavy- 
331]. The same is true for the 2 ^ and 2 operators. 


In the nonrelativistic basis, the I"*" and operators are distinguished by the presence 
of a Pauli matrix in the latter but not in the former. Terms in the NRQCD propagator 
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TABLE V: Spin splittings for IS*, 2S, IP, 2P and ID bottomonia. ^P represents the spin-averaged 


triplet P wave. Only statistical errors are shown for this worE 
“ Extracted from Eig. 5 of Meinel and Fig. 23 of HPQCD 
^ Extracted from Ref. 10]. 



this work Meinel [6, ^ 

HPQCD [9, JR, J^ experiment [30] 

T{lS)-r],{lS) 

57.97(16) 60.3(5.5)(5.0)(2.1) 62.8(6.7) 

62.3(3.2) [Belle=57.9(2.3)] 

T{2S)-r]h{2S) 

22(3) 

23.5(4.1)(2.1)(0.8) 26.5(1.6)(1.4) 

24(4) 

PP-Xbo{lP) 

35.7(5) 

36(3)“ 

40.0(1.2)“ 

40.4(5) 

l^P-Xbi{lP) 

8.3(4) 

8.5(1.8)“ 

7.4(6)“ 

7.1(4) 

l^P-h{lP) 

0.8(3) 

0.04(93)(20) 

2 .0(6)“ 

0 .6(1.0) 

Xb2{lP) - l^P 

12.1(3) 

12.4(1.4)“ 

12.4(5)“ 

12.3(2) 

2^P-Xbo{2P) 

24(9) 

- 

- 

27.7(7) 

2^P-Xbi{2P) 

6(5) 

- 

- 

4.8(5) 

23P - hb{2P) 

0(4) 

- 

- 

0.4(1.3) 

Xb2i2P) - 23P 

9(4) 

- 

- 

8.4(3) 

T2(1P) - T(1P) 

11(3) 

- 

17(6) 

- 

T3(1 P)-T2(1P) 11(2) 

- 

18(5) 

- 

T3(1P) - T(1P) 

22(3) 

- 

34(8) 

- 

r]b2 {lD)-T2 {lD) 1.6(1.5) 

- 

5(8)^ 

- 


that contain a Pauli matrix are also proportional to an odd power of the quark momentum. 
The cross correlator of a Pauli matrix operator with a non-Pauli matrix operator will be 
proportional to an odd power of momentum, and should be zero in the ensemble average 
because of spatial lattice symmetry. We verihed that cross correlators of the 1"' and I’*"''' 
operators for bottom mesons are statistically consistent with zero at all Euclidean times. 
This conhrms that the two operators are orthogonal, but it does not provide information on 
the mixing of the physical states within each operator. The same is true for the 2 ^ and 
2 operators. 

In practice the and operators appear to plateau at different energies, given a 
limited extent in Euclidean time. Even though they should contain the same ground state, 
they have a mixing with an excited state that is very close in energy to the ground state. 
They each have a different mixing with these two states, which can give the appearance of 
different ground states in a practical lattice study. These false ground-state signals would 
each be larger than the physical ground state and smaller than the physical excited state. 
A multiexponential £t is unable to distinguish the two physical states given the precision 
of our data. An application of the variational method to a correlator matrix can separate 
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these states (see, for example, Ref. [3^), but that is beyond the scope of this project. We 
will state the results from both I’*' operators (and both 2“ operators), using a prime for the 
heavier of the pair, and acknowledge that there is an unresolved mixing. 

For Bf. mesons, the S-wave correlator is fitted with five exponentials, while the P and D 
waves use four exponentials. The D-wave Tf operators have an overlap with the S-wave 
R*. The smallest energy is interpreted as the B* ground state and next smallest energy as 
the D-wave ground state. For Bg and R, the S- and P-wave correlators are htted with four 
exponentials. 

For bottom mesons, the additive NRQCD mass of Eq. flTH]) takes the form 


o f he 

me = - 


T-isim TTisim 

H/ - Il/'r 

2 ^ 


(23) 


and the lattice spacing from Eq. flT^ is used. Absolute masses are not calculated for Be 
mesons. The absolute mass for a Be meson contains large discretization effects because the 
mass of the charm quark is large compared to our lattice cutoff, so absolute masses will 
not be studied in this work. Instead, we focus on mass differences among Be states because 
these are expected to be close to their physical values [24 1. 

The spectrum of Rc mass differences with respect to the lightest Be state for IS, 2S, IP, 
2P and ID is shown in Fig. |6l The = 0“ radially excited S-wave state agrees with a 
recent observation by the ATLAS Collaboration 3^. First lattice results are given for the 
Rc D-wave ground states and the P-wave radial excitations. The Be results are all under or 
very near the BD threshold. Precise Rc spin splittings for IS, 2S, IP and IR are given in 
Table I VI I Note that spin splittings are only given between spin 0 and spin 2 for P waves, 
and spin 1 and spin 3 for D waves, because of the unresolved mixing between the two spin-1 
P waves and the two spin-2 D waves. 

The spectrum of Bg masses for IS, 2S and IR is shown in Fig. [71 Bg spin splittings 
are given in Table I VII The Bg radial and orbital excitations are either very near or above 
the threshold for breakup into BK or B*K. Therefore, the 2S and IP Bg states have the 
possibility to mix with these two-meson scattering states. A cross-correlation matrix of Bg 
with two-meson operators and application of the variational method, which was done in 
Ref. [l^, would be necessary to analyze the mixing and account for its effect. 

Since the = 2+ state does not decay to BK via S wave, and the S-wave decay to B*K 
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FIG. 6: Spectrum of meson mass differences with respect to the lightest B^ meson. The red 
band is an experimental value 35|]. Black points with errors bars are lattice data with statistical 
errors only. Grey bands are combined statistical and systematic lattice spacing uncertainty, added 
in quadrature. Numerical values are given in Table IVTIIl 


is suppressed for the = I"*" state with j = 


19 


, the two-meson decay must contain 


nonzero momentum. Therefore the hnite-volume thresholds for these two states are higher 
than the continuum threshold. The same thing applies for the S-wave radial excitations, 
which are forbidden to decay to by S wave. Given the minimum nonzero momentum 

for the lattice ensemble used in this study, the threshold is around 5.96 GeV for 5*2 and 
B*{2S), and 6.01 GeV for 5'^ and Bs{2S). Recall the systematic uncertainty due to lattice 
spacing: if the PAGS-GS lattice spacing from Eq. flTTll is used, then the B *2 — Bs and B'^^ — Bg 
mass differences agree with Ref. 0 and with experiment. 

The 0’*' and states for Bg have not been observed by experiment. Our results a gree 
with Ref. 16| which did not calculate mixing with two-meson states, but differ from Ref. 191 


which included two-meson states explicitly, as shown in Table I VI I This is another indication 
that including two-meson operators is a necessary step for obtaining reliable masses near 
and above threshold. 
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TABLE VI: Comparison of lattice results for B, Bg and Be- All quantities are in units of MeV/c^. 
Only statistical errors are shown for this work. 

“ Mixing with two-meson scattering states was not calculated. 

^ There is unresolved mixing between the two I"*' states. 



this work HPQCD 
[16,18] 

Lang et. al. ALPHA 

Ua BD 

experiment 

[30,^ 

Bs-B 

84(5) 

84(2) 

81.5(4.1)(1.2) 88.9(5.7)(2.3) 

87.35(23) 

B* - B 

41(3) 

50(3) 

47(7)(1) 41.7(4.2)(3.2)(0.3) 45.0(4) 

B: - Bg 

46.9(5) 

52(3) 

47.1(1.5)(0.7) 37.8(3.3)(5.8) 

48.7(2.2) 

B* - Pe 

57.5(3) 

54(3) 

- 

- 

B*{2S)-B{2S) 

60(25)“ 

- 

- 

- 

B*{2S) - Bg{2S) 48(9)“ 

- 

- 

- 

B*{2S) - B^{2S) 35(2) 

29(26) 

- 

- 

^2 ^0 

148(17)“ 

- 

- 

- 

^s2 ^sO 

104(7)“ 

- 

142(26) 

- 

^c2 -^cO 

68(3) 

- 

- 

- 


39(10) 

- 

- 

- 

B{2S) - B 

617(42)“ 

- 

791(73)“ 

- 

B*{2S)-B* 

636(39)“ 

- 

- 

- 

Bg{2S) - Bg 

594(14)“ 

- 

566(57)“ 

- 

b;(2S) - b: 

595(15)“ 

- 

- 

- 

P,(2^) - P, 

568(6) 

616(19)(1) 

- 

565(4)(5) 

B-J2S) - B- 

545(6) 

591(18)(1) 

- 

- 

b:o 

5770(6)“ 

5752(16)(5)(25)“ 5711(13)(19) - 

- 

Bgi 

5822(5)“’^ 

5806(15)(5)(25)“’^ 5750(17)(19) - 

- 

B'g, 

5844(6)“’^ 

- 

5831(9)(6) - 

5828.7(4) 

b:2 

5874(7)“ 

- 

5853(11)(6) - 

5839.96(20) 

B^o - Be 

436(2) 

429(13) 

- 

- 

Bel - B* 

419(2)^ 

410(13f 

- 

- 

B'ei - Bl 

431(2)^ 

- 

- 

- 


The spectrum of B masses for IS", 2S and IP is shown in Fig.lHl The LHCb Collaboration 
has evidence for two states P(5840) and P(5960) in the mass range consistent with B{2S) 
and B*(2S) 37l |. The CDF Collaboration also has evidence for a state P(5970) in this energy 
range 38(]. B spin splittings are given in Table IVIi The B radial and orbital excitations 


are all well above the continuum Pvr threshold. As discussed above for the case of P^, the 


two-meson S-wave decays are suppressed for B[ [36| and are forbidden for B*{2S) 


and B{2S). Given the minimum nonzero lattice momentum used in this study, the threshold 
for P 2 and B*{2S) is 5.76 GeV, and for B[ and B{2S) is 5.81 GeV. A thorough treatment of 
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FIG. 7: Mass spectrum of the Bg meson. Red bands are experimental values. Black points with 
errors bars are lattice data with statistical errors only. Grey bands are combined statistical and 
systematic lattice spacing uncertainty, added in quadrature. Numerical values can be found in 
Table US 


the B meson spectrum, like that done in Ref. 


39| for the D meson P waves, would require 


the inclusion of two-meson operators. 


VII. CONCLUSIONS 

In this work, the use of free-form smearing has allowed a lattice study of several bottom 
mesons and bottomonium states. The method permits the user to build a source operator 
with any desired shape, and does not require gauge hxing. New to the present work is a 
“minimal-path” method for free-form smearing that reduces the computational cost by 0{L) 
on xT lattices. Moreover, any number of unrelated smearing choices for a particular field, 
i.e. Eq. ([3]), can be built by reusing a single sum over the shortest link paths, i.e. Eq. ([5]). 

Beginning with local operators designed to have the appropriate quantum numbers 
for each hadron, free-form smearing was used to produce a particular hydrogen-like wave- 
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FIG. 8; Mass spectrum of the B meson. Red bands are experimental values. Black points with 
errors bars are lattice data with statistical errors only. Grey bands are combined statistical and 
systematic lattice spacing uncertainty, added in quadrature. Numerical values are given in Table 1X1 


function shape in each case, with the radius and node positions tuned for optimal results. For 
each meson, the free-form smearing was applied to the bottom quark at the source, leaving 
the antiquark {u/d, s, c or b) unsmeared. Bottom quarks were handled with NRQCD and 
all other quarks were relativistic. 

A nearly complete picture of the bottom meson and bottomonium spectrum below the 
two-meson breakup threshold is given in Figs. [5l El [7] and [HI with hrst lattice results for 
bottomonium D-wave radial excitations. Be D-wave ground states and P-wave radial 
excitations. The results for states that are well below the breakup threshold (e.g., BB^*'> 
for bottomonium, B^*^D for Be) agree reasonably well with available experimental data and 
several predictions were made. 

Some states that were near or even above threshold are also calculated. These include 
the B and Bg S-wave radial excitations and P-wave ground states, the Be P-wave radial 
excitations and the bottomonium G-wave ground states. Where mixing with meson scat- 
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tering states may be important, a more complete analysis will be required to obtain robust 
predictions. This work, aimed at constructing correlation functions with good radial and 
orbital excitation signals, is a step in that direction. 

The original free-form smearing jd] was considered unusable for sink smearing. The 
speedup obtained through our minimal-path implementation gets us closer to the possibility 
of sink smearing and is an avenue for further development which would make the technique 
more generally applicable. There are also possibilities to combine this with other ideas like 
the variational method. A specific application for further development is the calculation of 
matrix elements, such as in three-point functions where typically one uses a fixed source-sink 
separation. 
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Appendix A: Additional tables 

Numerical values from Figs. O |6l [7] and [8] are provided in Tables I VI II IVIIIl HXl and 1X1 
respectively. 
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TABLE VII: Bottomonium simulation energies in lattice units and masses in units of MeV/c^. 
Lattice values in physical units for the J = 2, 3 states are dimensionally averaged over the different 
lattice irreducible representations. Only statistical errors are shown for this work. 


particle 


simulation energy 

mass [MeV/c^ 

this work experiment [29, .3d] 

Vb{lS) 

0 -+ 

0.23401(7) 

9402.3(3) 

9398.0(3.2) 

T(l^) 

1 -- 

0.25963(10) 

9460.30(26) 9460.30(26) 

Vb{2S) 

0 -+ 

0.497(3) 

9998(6) 

9999(4) 

T{2S) 

1 -- 

0.507(3) 

10020(7) 

10023.3(3) 

r]b{3S) 

0 -+ 

0.637(19) 

10314(44) 

- 

T(3^) 

1 ” 

0.646(23) 

10334(52) 

10355.2(5) 

h{lP) 

1 +- 

0.4540(5) 

9900.2(9) 

9899.3(1.0) 

Xboi^P) 

0 ++ 

0.4386(4) 

9865.3(1.0) 

9859.4(5) 

Xbii^P) 

1 ++ 

0.4507(5) 

9892.7(9) 

9892.8(4) 

Xb2{lP) 

2 ++ 

0.4599(6)e 

0.4596(6)t2 

9913.1(1.0) 

9912.2(4) 

hb{2P) 

1 +- 

0.595(5) 

10219(11) 

10259.8(1.2) 

Xbo{2P) 

0 ++ 

0.584(5) 

10194(12) 

10232.5(6) 

Xbi{2P) 

1 ++ 

0.592(5) 

10212 (11) 

10255.5(5) 

Xb2{2P) 

2 ++ 

0.600(6)i5 

0.598(6)t2 

10227(13) 

10268.7(5) 

Vb2i^D) 

2 -+ 

0.571(3)^5 

0.5693(15)t2 

10163(3) 

- 

T{1D) 

1 — 

0.5646(9) 

10150(2) 

- 

T^ilD) 

2 — 

0.5689(15) 

0.5697(17) 

10161(3) 

10163.7(1.4) 

TsilD) 

3— 

0.5728(18)^2 

0.5761(10)ri 

0.5730(16)r2 

10172(3) 


Vb2i2D) 

2 -+ 

0.712(11)£; 

0.693(18)t2 

10458(26) 

- 

T{2D) 

1 — 

0.675(17) 

10401(39) 

- 

T2(2D) 

2 — 

0.690(16)i? 

0.699(10)t2 

10447(25) 

- 

T,{2D) 

3— 

0.703(24)^2 

0.704(22)ti 

0.698(15)t2 

10459(37) 


hsilP) 

3+- 

0.654(4)^, 

0.655(4)^2 

10355(7) 

- 

Xbsi^F) 

3++ 

0.653(5) 

10350(12) 

- 

Vbii^G) 

4-+ 

0.749(4) 

10568(9) 

- 

T4(1G) 

4— 

0.760(9) 

10592(21) 

- 
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TABLE VIII: B^. meson simulation energies in lattice units and mass differences with respect to 
the lightest Be state in units of MeV/c^. Lattice values in physical units for the J = 2,3 states are 
dimensionally averaged over the different lattice irreducible representations. Only statistical errors 


are shown for this work. 


particle simulation energy 


this work experiment [30, 


35 


Be{lS) 0- 

0.73373(14) 

0 

b;(is) 1- 

0.75914(14) 

57.5(3) 

Bc{2S) 0- 

0.985(2) 

568(6) 

Bt(2S) 1- 

1.000(3) 

603(6) 

B:„(.1P) 0+ 

0.9266(8) 

436(2) 

BciilP) 1+ 

0.9445(9) 

477(2) 

B'elilP) 1 + 

0.9496(10) 

489(2) 

B’JIP) 2+ 

0.9564(13)^5 

505(3) 


0.9569(13)t2 


b;„(2p) 0+ 

1.114(9) 

861(22) 

Bci{2P) 1+ 

1.121(10) 

875(24) 

B;,(2F) 1+ 

1.120(11) 

875(26) 

B*A^P) 2+ 

1.127(12)^5 

881(27) 


1.120(13)r2 


b;{id) 1- 

1.0674(18) 

755(5) 

Bc2{1D) 2- 

1.074(3)e 

775(9) 


1.077(5)t, 


B'JID) 2- 

1.072(14)^ 

777(16) 


1.080(7)t2 


b;jid) 3- 

1.089(8)^, 

794(8) 


1.085(2)ti 



1.083(6)t2 



0 

565(6) 
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TABLE IX: Bg meson simulation energies in lattice units and masses in units of MeV/. Lattice 
values in physical units for the J = 2 states are dimensionally averaged over the different lattice 
irreducible representations. Only statistical errors are shown for this work. 


particle simulation 

energy mass 

this work 

MeV/c^; 
experiment [30j 

Bg{lS) 0- 0.4130(3) 

5370.9(1.6) 5366.77(24) 

Bl{lS) 1- 0.4337(4) 

5417.8(1.8) 5415.4(2.3) 

Bg{2S) 0- 0.675(6) 

5965(14) 

- 

Bl{2S) 1- 0.697(7) 

6013(15) 

- 

B*q{IP) 0+ 0.590(2) 

5770(6) 

- 

Bsi{lP) 1+ 0.612(2) 

5822(5) 

- 

fi'i(lP) 1+ 0.622(3) 

5844(6) 

5828.7(4) 

B:^{IP) 2 + 0 . 635 ( 3 )e 

5874(7) 

5840.0(2) 

0 . 636 ( 3 ) t 2 




TABLE X: B meson simulation energies in lattice units and masses in units of MeV/c^. Lattice 
values in physical units for the J = 2 states are dimensionally averaged over the different lattice 
irreducible representations. Only statistical errors are shown for this work. 


particle simulation 

energy 

mass [MeV/c^ 

this work experiment [30, ^ 

P(l^) 0“ 0.3757(20) 

5287(5) 

5276.26(17) 

B*{IS) 1- 0.3937(19) 

5327(4) 

5325.2(4) 

B{2S) 0- 0.648(20) 

5904(44) 

- 

B*{2S) 1- 0.675(18) 

5964(41) 

- 

B*{IP) 0+ 0.544(4) 

5667(10) 

- 

Pi(lP) 1+ 0.575(3) 

5738(8) 

- 

P((1P) 1+ 0.590(4) 

5771(10) 

5727.7(1.6) 

P*(1P) 2+ 0.609(7)^ 

5815(14) 

5739.4(5) 

0.610(6)r2 
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